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ABSTRACT 

Context. Hydrodynaniical model of the ^-mechanism in a purely radiative case. 

Aims. First, to determine the physical conditions propitious to K-mechanism in a layer with a configurable conductivity 

hollow and second, to perform the (nonlinear) direct numerical simulations (DNS) from the most favourable setups. 

Methods. A linear stability analysis applied to radial modes using a spectral solver and DNS thanks to a high-order 

finite difference code are compared. 

Results. Changing the hollow properties (location and shape) lead to well-defined instability strips. For a given position 

in the layer, the amplitude and width of the hollow appear to be the key parameters to get unstable modes driven 

by K-mechanism. The DNS achieved from these more auspicious configurations confirm the growth rates as well as 

structures of linearly unstable modes. The nonlinear saturation follows through intricate couplings between the excited 

fundamental mode and higher damped overtones. 
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1. Introduction 

Since the beginning of the studies concerning Cepheids, 
it is well known that convection occurs in Cepheids' en- 
velopes, and thus changes pulsation propertie s (e.g. the re- 
views of lGautschv fc Saiolll996l : lBuchleilll99l . The coldest 
ones, which are located next to the red edge of the instabil- 
ity strip, have the more extended surface convective zones. 
However, during many years Cepheids' oscillations models 
have used the so-called "frozen-in convection" approxima- 
tion whic h claims that convective flu x perturbations are 
negligible (jBaker fc Kippenhahnl 11963 ). Such kind of mod- 
els well predict the blue edge of the instability strip but 
fail to explain the red edge as in this case, the strong exist- 
ing couplings between the convection and the oscillations 
are not taken into account. This discrepancy becomes obvi- 
ous with the accumulation of accurate observations which 
show a narrower instability strip than the theoretical one, 
i.e. modes which are linearly unstable in the models are 
located outside the observational strip (JYecko et al.lll998l 
hereafter YKB98). 

The main theoretical difficulty comes from the fact that 
convection plays a crucial role on the pulsations while we 
know that convection i tself remains roughly described b y 
mixing-length theories (|Vitensdll953tlBohm-Vitenselll958f ). 
However, several attempts in the direction of time depen- 
dent convection models (TDC) have been developed (e.g. 
[U^o 1967: Gough 1977HStellingwerjll982[ ). Recent studies 
(YKB98; iBono et al.lll999D re ly on SteUingwerf's convec- 
tion model ( Stellingweril Il982f ) or very similar newer ap- 
proaches to compute linear and nonlinear time evolution of 



ampl i tudes of modes (iKuhfufil 1986t iGehmevr fc Winkleij 
119921: IWuchterl fc FeuchtingeJll998[ ). The major problem 
of TDC is the choice of the many free parameters intro- 
duced by the convection modeQ- As these parameters are 
not theoretically well determined, one should adjust them 
by fitting the observations. 

Another way to study the convection-pulsation interac- 
tion is to achieve (nonlinear) direct numerical simulations 
(DNS) of the whole hydrodynamical problem. The final aim 
of our work is to realise such kind of simulations in 2-D 
and 3-D where a convective zone will be coupled with a 
radiative one and unstable radial acoustic modes will be 
self-consistently excited by K-mechanism. However, as DNS 
arc highly time consuming, it is necessary to get in a first 
step the appropriate initial conditions. That is why we have 
tried to determine in this paper the physical conditions pro- 
pitious to an excitation based on K-mechanism. 

jEddingtonI (|1917[ ). and then IZhevak"h3 (|1953[ ) and[Co3 



()1958D have introduced a mechanism linked to the opacity 
in ionisation regions, the K-mechanism, where k denotes 
opacity (see also the review of Zhevakin 1963,). They have 
shown that Cepheids' radial oscillations are driven by a 
thermal heat engine as radial pulsations have to be main- 
tained thanks to a sustained physical process. That is why, 
they imagined a Carnot-like thermodynamic cycle which 
stores heat during compression phases while releasing it 
during the decompression ones. This mechanism, now called 
the Eddington's valve, can only occur in specific regions of 
a star where the opacity varies so as to block the radia- 
tive flux during compression phases. Now, opacity tables 
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^ e.g. the seven coefficien ts [qc, Qt, cty , ah, QJs, ad, "p] in 
YKB98 or the eight ones in lKoUath et all lj2002f ). 
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such as the OPAL ones show strong increases in ionisation 
regions of main (i.e. Hydrogen and Hehum) or heavy ele- 
ments that are common ly called "ionisation bumps" (e.g. 
ISeaton fc BadneUl [2Q0J) . As a consequence these ionisa- 
tion zones are locally responsible for modes amplification. 
Moreover, beyond this criterion on opacity, these ionisation 
regions have to be located in a very precise region of a star. 

Indeed, if they are too deep or too close to the surface, 
the driving they cause can be balanced by the damping oc- 
curring in other regions. Therefore, an efficient ionisation 
region has to be located in a specific place, called the tran- 
sition region, which marks the shift from a quasi-adiabatic 
interior to a strongly non-adiabatic surface. In classical 
Cepheids oscillating on the radial fundamental mode, the 
overlap of this transition region with the ionisation one is 
around 40 000 K, which corre sponds to the temperature o f 
the helium second ionisation ([Baker fc KippenhahnlflQSa l. 
For first overtone Cepheids, things are more intricate as 
one must take into account first, the respective position 
between this Hell ionisation region and the nodal line and 
seco nd, the Hel/H reg ion which also contributes to the driv- 
ing (|Bono et al.lll999l) . 

However, the location of these opacity bumps are not 
solely responsible of the acoustic mode destabilisation. A 
careful treatment of the K-mechanism in Cepheid stars 
would also involve the possible dynamical couplings with 
convective zones or, say, metallicity effects through a realis- 
tic equation of state and opacity tables. The corresponding 
physics is actually not fully understood from a theoretical 
point of view and the purpose of our work is to sufficiently 
simplify the hydrodynamical approach while keeping at the 
same time the leading order phenomenon responsible of the 
instability, that is, the opacity bumps location. 

That's for why we adopt a fully radiative layer of a per- 
fect gas in which a ionisation region is represented by a "hol- 
low" in radiative conductivity, corresponding to a "bump" 
in opacity. Strictly speaking, the layer stability will depend 
on both temperature and density variations, as the radia- 
tive conductivity is a function of these two physical quan- 
tities. However, as the opacity strongly depends on tem- 
perature, this state variable mainly controls the instability. 
As our final aim is to realise an hydrodynamical study of 
K-mechanism, we will therefore neglect any dependence of 
radiative conductivity on density and thus, our conductiv- 
ity profile is merely a function of temperature. The inferred 
advantage is to allow easy changes in the different param- 
eters of the ionisation region by setting both position and 
shape of the conductivity hollow (i.e. its slope, width and 
amplitude). As a consequence, it is possible to investigate 
a complete parametric study of K-mechanism in order to 
determine precisely the physical conditions required by the 
instability. 

We first consider the radiative and hydrostatic equilib- 
ria of our layer with an appropriate conductivity profile. 
By adjusting both the density value at the top of the layer 
and the fiux at its bottom, we then obtain a transition re- 
gion in the middle of the computational domain. Secondly, 
we investigate the linear stability analysis by solving per- 
turbation equations for the oscillations. We thus obtain the 
whole spectrum, and can sort out unstable modes from sta- 
ble ones. Therefore, we are able to check the relevance of 
the transition region concept by making every parameters 
of the conductivity hollow vary. 



The first main result of our linear stability study is the 
confirmation of the underlying conditions defining the tran- 
sition region. With different parametric studies, we obtain 
instability strips corresponding to the fundamental mode, 
as both a minimum hollow width and amplitude are needed 
to obtain unstable modes. These results are interpreted 
thanks to the work integral to exactly determine the loca- 
tion of the driving zone. As a consequence, this parametric 
study of K-mechanism provides us the physical quantities 
responsible for the instability. 

Secondly, these auspicious conditions needed to drive 
the fundamental mode constitute the starting point of 2- 
D DNS. Indeed, we check the growth rates as well as the 
structure of the linearly unstable modes by performing di- 
rect simulations until the nonlinear saturation of modes. 

In f}2l we first introduce the general oscillation equations 
while the different conditions leading to the instability are 
determined in ^ In 21 we present our hydrodynamical 
model. Linear stability analysis and DNS results are thus 
given and compared in fJ5]and fJSl respectively, before con- 
cluding in 121 

2. The general oscillations equations 

Our system is composed by a 2-D plane parallel layer of 
width d. z denotes a cartesian coordinate, pointing upward 
on the contrary of the constant gravity field g. The gas is 
assumed to be monatomic and perfect, so its equation of 
state is given by 



•p — R^pT and 7 ; 



(1) 



where p, p and T respectively denote pressure, density and 
temperature, i?* is the ideal gas constant and 7 the ratio 
of specific heats Cp and c^ . 
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Fig. 1. Scheme corresponding to our model. Gravity is 
pointing downward, contrary to the vertical coordinate. 
The (blue) curve represents the radiative conductivity pro- 
file that we are going to discuss further whereas large (red) 
arrow expresses the radiative flux entering in the bottom 
of the layer. 

We are interested in small perturbations around an equi- 
librium state, that is, any physical quantity is expanded 
around its mean value Xo(r) as 

X{r,t)^XQ{r)+X'{r,t) with X'/ATo < 1, (2) 
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where X' is an Eulcrian perturbation. Linearized continu- 
ity, momcntuni and energy equations for the perturbations 
in a non-adiabatic case are given by (e.g. lUnno et al.lll989[ ) 



f Ap' 



XT' 



— po div M — It • V In po 

Po Pa 

+j/ [ Ait+ -VdivM + 2S'- Vhipo 



■ divF' - (7 - l)To divM - u ■ VTq, 



PqCv 



(3) 

where u, p' T' and F' denote the velocity, density, tem- 
perature and radiative flux perturbations, respectively. The 
kinematic viscosity v is supposed to be constant and S de- 
notes the traceless ratc-of-strain tensor given by 



5,, 



duj 
dx-i 



9uj 
dxi 



-Sij divu 



(4) 



We seek normal modes with a time-dependence of the form 
exp(Ai), where A = r+iuj. The real part r denotes the grow- 
ing (r > 0) or damping (r < 0) rate whereas the imaginary 
part tu denotes the frequency. We now assume that the layer 
is fully radiative which, under the diffusion approximation, 
leads to the following expression of the radiative flux per- 
turbation 



F' = -/voVT'-/\'VTo, 



(5) 



where Kq denotes the radiative conductivity and K' its 
Eulerian perturbation. This perturbation K' can be related 
in a general way to the temperature one by 



= K,T 

Ko To 



ith K-T = 



d In Kr 



cflnTo 
Finally, we impose the following boundary conditions 



(6) 



U2 = for z 
dT' 



[0, d]. 



(7) 



. dz 



for z == and T' = for z 



which correspond to rigid walls at both limits of the do- 
main, a perfect conductor at the bottom and a perfect in- 
sulator at the top. 

3. Conditions for instability 

3.1. A first condition derived from the work integral 

We recall that the main aim of this work is to clarify 
the favourable conditions which may sustain unstable ra- 
dial modes in a plane parallel layer. An advisable mean to 
study the physics of such instability is the work integral. In 
the following, we are going to assume that transformations 
are quasi-adiabatic. This approximation is sufficient in the 
deeper layers of a star but becomes no longer valid near its 
surface where non-adiabatic effects dominate. Nevertheless, 
this approximation is useful when one wants to have an idea 
of the stability properties of an oscillation mode. 



Using the work integral formalism in the quasi-adiabatic 
limit, we demonstrate in Appendix lA.ll the following ex- 
pression for the damping or growing rate of an eigenmode 



5R 



(^_l)^divF'dV" 
V P 



(8) 



H^PodV 



where the symbol 3? means the real part. We thus see that 
the sign of r only depends on the numerator which, under 
the diffusion approximation, Eq. ([5]), leads to 



r ex 3fi 



Po 



diY{K'VTo + KaVT')dz 



(9) 



The main driving term in this expression is K''^Tq because 
it represents the dynamical variation of opacity during an 
oscillation cycle, which is the cause of K-mechanism. We 
thus neglect Kq'VT' to get 



3? 



^ div(i^'VTo)dz 
Po 



(10) 



As we can see in Fig. [5]d, the equilibrium temperature is 
an almost linear function of z, except in the vicinity of 
the conductivity hollow. Therefore, VTq is almost constant 
thus we have 



T (X 3ff 



— VToVK'dz 
Pq 



(11) 



With Eq. I®, we have K' = Kq{Tq)1CtT' /Tq. As we are 
interested in low order modes we assume that V-R"' is dom- 
inated by Kq{To)T' /Tq'V ICt ■ The same approximation will 
be done in Eq. (|17p . We then obtain 



r ex 3fi 



^ VToA'o(To) — V/Crdz 



Po 



Tn 



(12) 



Let us now consider a compression phase corresponding 
to Sp*/po > and T'/Tq > 0. As VTq < 0, a necessary 
condition to obtain unstable modes with r > is then 



dK-T 
dz 



<0. 



(13) 



In variable stars, this condition can occur in ionisation re- 
gions. Indeed, opacity tables clearly show that these re- 
gions are ass ociated with large "bumps" in opaci ty (e.g. 
the review of ICarsonlfTOTJISeaton fc Badnel]||2004[ ). With 
Eq. (|13[) . one can obtain unstable modes if driving pre- 
vails over damping. This result has the following physical 
meaning: if the radiative conductivity decreases during a 
compression phase, then some part of the radiative flux is 
blocked and some energy is stored during compression con- 
tributing to increase the ratio of ionised matter. During 
the decompression phase, this extra energy is released un- 
der mechanical work to the environment and can excite the 
oscillations. It works as Eddington predicted when he imag- 
i ned a thermal he at engine to justify Cepheids' oscillations 
(|Eddingtonlll917t) . 
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3.2. A second condition derived from the so-called 
"transition region" 

In the previous paragraph, we have recalled a necessary 
condition to obtain unstable modes (Eq. [13]). Nevertheless, 
the demonstration was made under the quasi-adiabatic ap- 
proximation which fails near the surface. As a consequence, 
we cannot know if the driving caused by ionisation regions 
can prevail over other damping regions. Indeed, ionisation 
regions where Eq. (J13p is satisfied are very thin compared to 
the whole atmosphere of a star and thus the influence of the 
damping regions on the instability remains questioning at 
this stage. He reafter, we will essentially follow Cox's d emon- 
stration fe.g. ICox|[T980l : IChristensen-Dalsgaardll2003[) . 

In Lagrangian variables, the energy equation can be 
written as 



— If the considered point is too deep in the star, the local 
thermal timescale is longer than the dynamical one and 
5* 3> 1. As a consequence, pulsation has no influence on 
the energy release and the adiabatic approximatioro is 
well-suited in this case. 

— On the contrary, if the considered point is next to the 
surface, then Am is very small and ^P ^ 1. As a conse- 
quence, the energy content (and the corresponding mass 
content) is too weak to influence the radiative flux and 
the luminosity perturbation is said to be frozen in. It 
means that the radiative flux perturbation doesn't vary 
anymore in this region. 

— Between these two regions, '^ can be 0(1) which defines 
a zone called the transition region separating the adi- 
abatic interior from the strongly non-adiabatic surface, 
that is, 



d 
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where the symbol 6 means Lagrangian perturbations. As 
we consider a 1-D box, we simply have 
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(15) 



Let us integrate between a given height z and the surface 
to get 



A 



F' 

'F~ 



bot 



n. 



{CyTp) 

Fhot 



podz' d 



ST 



(7-1) 



S_p 
Po. 



(16) 

where A(F'/Fbot) is the variation of -FV-Fbot between the 
considered point and the surface, {c^Tq) is an average of 
CyTp over this region and 11 is a characteristic dynamical 
timescale (the pulsation period of the fundamental mode, 
typically). As we are principally interested in the study of 
low order modes (i.e. the fundamental one) , we assum e that 
the eigenvectors are almost constant (e.g. ICoxlll980f ). As a 
consequence, one obtains 



5T , ^.Sp 
— -(7-1) — 
J-o Pa. 



ST , ^.Sp 

— -(7-1) — 
Jo Po 



(17) 



where {.)z expresses an average of the eigenf unctions. 
Noting that 



{cyTo)Am ^ 
UL ^ ' 



(21) 



iKing fc Coxl (jl968f ) have first determined the temperature 
associated with this transition region for different Cepheids' 
models. For the fundamental mode, they have obtained 
Ttr — 40 000 K while the transition lies nearer the sur- 
face for higher order modes. 

The corresponding positions of ionisation regions and 
transition ones are therefore crucial for the instability 
(jGiUiland et al.l [l998f ). Indeed, if they overlap -it corre- 
sponds to the instability strip- the bottom of the ionisa- 
tion region strongly contributes to driving because it acts 
in a quasi-adiabatic place. On the contrary, its top is in a 
strongly non-adiabatic region where the luminosity is frozen 
in. As a consequence, the zone over ionisation region is not 
damping and dri ving prevail s in this case: modes become 
unstable (see e.g. ICoxl (|1980( ) for a more detailed descrip- 
tion). 

We then sum up the conditions propitious to instability 
as 

< and Ttr ~ Tionisation- (22) 



dz 

With the conditions given in Eq. ([H]), one thus obtains 
^ionisation — 40 000 K. This temperature corresponds to 
the second Helium ionisation which is known to be m ainly 
responsible for the driving of modes in Cepheids (see ICoxl 
:196a . 



leads to 



where 
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podz' = — ^ and L 



f F' 



V-Fbo 
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^ 



ST Sp 

77r-(7-l) — 
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(18) 4. Our hydrodynamical model 

4.1. The choice of the radiative conductivity profile 



(19) 



(20) 



^ has the following physical meaning: it represents the ra- 
tio of the thermal energy embedded between the considered 
point and the surface to the energy radiated during an os- 
cillation pcrioqj. A second instability criterion can then be 
derived from the value of this quantity ^P: 

■^ 't can also be interpreted as the ratio of the local thermal 
timescale to the dynamical one. 



Strictly speaking, radiative conductivity Kq depends both 
on tempe rature and density as diffusion approximation 
gives (e.g. iMihalas fc Weibel MihaTaill984t ) 



Ko 



16ctT3 
inp 



(23) 



where a denotes Stefan-Boltzmann's constant. Kramer's 
laws constitute good approximations of opacity laws with 



One can notice that the "adiabatic" approximation means 
here that deeply in the layer the mode has an adiabatic be- 
haviour whereas next to the surface it remains strongly non- 
adiabatic. 
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K oc p"r- 



(24) 



where, for ex ample, n = 1 and s = 3.5 for the free-free 
opacity (e.g. [Hansen fc Kawalen 11994 ). As shown in the 
Introduction and in ^ we simply consider here a con- 
stant radiative conductivity on which a hollow representing 
a bump in opacity is added. This simple approach is an easy 
way to reproduce the physical conditions propitious to n- 
mechanism while keeping the ability to quickly change the 
hollow amplitude, slope or width and achieving so an effi- 
cient parametric study of the instability. As a consequence, 
the radiative conductivity is given bjQ 



A'o(T) = K^ 
with 



l + A 



-7r/2-harctan(crT+T~) 



A^ 



K — K 



7r/2 -f arctan(cre^) 



-^ — -^ -'bump ^ ^5 



(25) 



and 



~ ^^bump is the position of the hollow in temperature. 

— i^niax and iCniin ^rc the conductivity extrema, A being 
the corresponding relative amplitude. 

— a represents the slope of the hollow. 

— e represents the half of the full width at half maximum 
(i.e. FWHM/2) of the hollow. 

Examples of conunon values of these parameters are given 
in Fig. [21 
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Fig. 2. Influence of the hollow parameters on the conduc- 



tivity profile for K^ 



10-2 and Tb 



3.5: amplitude 



A (a) , width e (b) and slope a (c) (here the abscissa denotes 
a dimensionless temperature based on the surface one, i.e. 

J- / J- surface/ 



4.2. The equilibrium setup 

Hydrostatic and radiative equilibria in the diffusion limit 
are given by 

r Vpo = Pag 

(27) 
idiv[A-o(To)Vro]=0. 

Assuming a constant radiative flux i^bot at the bottom of 
the layer leads to 



dz 

dJh 
dz 



-pog 



Fhot 



(28) 



Let us choose the depth of the layer d as the length scale 
and the temperature at the top of the layer Ttop as the tem- 
perature scale, i.e [z] = d and [T] = Ttop- Moreover, the top 
density is chosen as the density scale {[p] = ptop), velocity 
(26) is given in units of -^CpTtop, gravity in units of CpTtop/d, 
pressure in units of ptopl^^], while radiative conductivity 
is given in units of ptopCpd[u\. Corresponding radiative flux 
)'^/2. The dimensionless equations 



unit is then ptop{cpTtop 
then become 



dlnpo 
dz 
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The set of equations ()29p can be written in matrix form as 

Alp = B{ip), (30) 

where ijj = {pq^Tq)^ \s the equilibrium field vector and A, B 
are differential operators. One notes that the RHS, B{xp), 
depends on the eigenvector itself through the terms I/Tq 
and l/Ko(To), which means solving a nonZmear problem. 

Finally, tildes on equilibrium fields emphasise dimen- 
sionless quantities but they will be dropped for clarity in 
the following. 

4.3. About Schwartzscliild's criterion in our model 

It is important to keep our box entirely radiative. That is 
why, we must take c are of Schwartzschild's criterion (e.g. 
IChandrasekhaHll96lh 



dTo 
dz 


< 


dTo 
dz 



with 



adi; 



dTn 



dz 



adii 



9_ 

Cr) 



(31) 



If this inequality is ensured, then the layer is fully radiative. 
In our dimensionless units, this condition becomes 



'^ This arctan-profile may appear quite intricate, compared e.g. 
to a gaussian-like one, but it allows us to change almost indepen- 
dently the hollow parameters (i.e. changing its amplitude while 
keeping its width). 



dTn 



dz 



<9, 



(32) 



meaning that large variations in temperature through the 
domain necessarily require large values for the dimension- 
less gravity field g. 
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4.4. Back to the radial oscillations equations 

To simplify the general system of oscillations equations ^ 
we first define the following new variables 



R = p'/po and e = T'/Ta. 



(33) 



Then, by using Eq. (|29|) and substituting the conductivity 
equation into the energy one (see Appendix IBJ) . one gets 



XR = — div u — u ■ 'Vlnpo 

Am = - — V{e + R)-0g 
Po 



-J/ ( Ait+ -Vdivit + 2S'- Vlnpo 



\e = — A,(i^oTo0)-(7-l)divM-M-Vlnro, 

PoCvJ^O 

We then adopt the same dimensionless quantities used in 
the equilibrium equations (P5|) 



4.5. The numerical methods 

We solve the two linear algebra problems ((30)l and 
(I371I38I) using the L SB code (Linear Solver Builder, 
IValdettaro et al.l l2007t ) . Both problems are discretized on 
the Gauss-Lobatto grid associated with Chebyshev's poly- 
nomials leading to two distinct numerical problems: 

1. the equilibrium model: the computation of the equi- 
librium structure requires to solve a nonlinear prob- 
lem. One way to do that is to use the so-called 
"Picard's method" based on the fix ed-point algorithm 
(jHairer et al.lll993l : [Fukushimalll997f) . It consists in solv- 
ing our set of first-order ordinary differential equations 
by successive iterations, that is, we advance 



Atp„+i = B{ip„). 



(39) 



This scheme converges quite well provided that the ini- 
tial guess ■00 is not "too far" from the solution and that 
nonlinearities are weak. 
2. the eigenmodes: we can compute the whole spectrum 
of complex eigenvalue s A using the QZ algorithm 
(jMoler fc Stewartlll973[ ) or just compute a pair (A, xp) 
around a given guess of A using the iterative A rnoldi- 
Chebyshev algorithm (IArnoldilll95it [Saadlfl99l) . 
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d In To 
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(35) 

where u^ denotes the vertical velocity, xo = Kq/ poCp the 
radiative diffusivity and the viscous dissipative term 2?^ is 
given by 



V,, 



4 f duz dlnpoduz 
3 \ dz^ dz dz 



(36) 



This system ()35|) may formally be written as a generalised 
eigenvalue problem 



Alp = XBip, 



(37) 



where A = r + icj is the complex eigenvalue associated with 
the eigenvector ip = {R,Uz,9)^ while A, B denote differ- 
ential operators. 

Finally the set of boundary conditions ([7]) written for 
(uz, 9) is 



( u,^0 for z = [0, 1] 
d9 din To 



dz dz 

9 = for z 



9 = for z 
1. 



(38) 



5. Results 

In Eq. (|32p . we have shown that the temperature contrast 
across the layer (and also the associated density and pres- 
sure ones) is limited by the gravity field g. It may how- 
ever be judicious not to restrict ourselves to small density 
contrasts as the mass involved between the conductivity 
hollow and the surface, that is Am, enters in one of the 
two favourable criteria for the instability, see Eq. (HI]). We 
will therefore consider in the following a "convenient" value 
5 = 7. 

Moreover, we have to keep as small as possible the ra- 
diative conductivity as we want to avoid too large values 
of diffusivity xoj that is, we do not want to deviate too 
much from conditions of the applying of quasi-adiabatic 
developments. We therefore choose Kj^^x = 10^^. As a con- 
sequence, the system (j29|) combined with Schwartzschild's 
criterion ([5^ restrict the possible values of the imposed 
bottom flux and we thus set Tbot = 2 x 10^^ in the follow- 
ing. Using mild spatial resolutions, we choose a conservative 
value for the kinematic viscosity, that is, z^ = 5 x 10^^. 

5.1. Computation of equilibrium fields 

To compute the equilibrium fields from second-order sys- 
tem (j29p . we must choose two different boundary condi- 
tions, one on temperature and one on pressure. Without 
loss of generality, we first set To = 1 at the top of the layer. 
Second, as we are interested in having the transition re- 
gion of the fundamental mode located in the layer middle, 
Eq. (PT]) has to be satisfied at that place. We have already 
said that Tbot = 2 x 10^^ and A'max = 10^^. Temperature 
value is fixed by Tbot, AT^ax and its boundary condition. 
The pulsation period of the fundamental mode is roughly 
expressed by 



U 



2d 



(40) 
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Fig. 3. a) Three different conductivity hollows: Tbump = 
2.8 (dashed red line), Tbump ~ 2.1 (solid green line) and 
Tbump = 1-7 (dotted blue hue), with A = 70%, e = 0.4 and 
cr = 7; b) corresponding temperature profiles; c) density 
profiles; d) equilibrium profiles for dK,T/dz. 



where (cs), denoting the mean sound speed, is related to 
the temperature contrast across the layer. Thus, the only 
parameter that we can change in Eq. (|21[) is Am. This value 
is directly linked to the boundary condition on density. In 
order to have a suitable ATn-value in the box middle, we 
therefore decide to take po = 2.5 x 10"'^ at the top of the 
layer (corresponding to InPo = —6.9). 

Once chosen the values of the hollow parameters - 
Thn-mp,-A, e and a- we then solve the problem ([29]) on the 
Gauss-Lobatto grid to obtain the three equilibrium fields 
To, Po and po- With these fields, it is possible to compute 
any physical quantities needed in the oscillations equations 
([55]). e.g. JCt or the diffusivity xo- Some examples of ob- 
tained equilibrium fields are given in Fig. [3] for three differ- 
ent positions of the conductivity hollow. 

5.2. Computation of eigenmodes 

With the different equilibrium quantities, we solve the 
eigenvalue problem p5p to get the whole spectrum of eigen- 
values A. Among them, it is possible to obtain the nearest 
eigenvalue from a guessed one using the iterative Arnoldi- 
Chebyshev algorithm. As an example. Fig. 3] brings out the 
eigenfunctions associated with the unstable fundamental 
mode. 

In order to check the convergence of this mode, we also 
compute the spectra corresponding to the different eigen- 
functions, see Fig. [S] 

Finally, in order to determine the influence of the con- 
ductivity hollow on that fundamental eigenvector, we com- 
pute three eigenmodes with the different values of Tbump 
previously used in Fig. [3l The result is displayed in Fig. [6] 
for the temperature perturbation 9, where the case of a 
constant radiative conductivity profile (i.e. without a con- 
ductivity hollow) is superimposed: 
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Fig. 4. Eigenfunctions (i?, u^, 6) for the unstable fun- 
damental mode, corresponding to a hollow with Tbump = 
2.1, A = 70%, e = 0.4 and a = 7. The equilibrium setup 
used to compute this mode is the one displayed by a solid 
green line in Fig. [31 
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Fig. 5. Chebyshev's spectra of the different fields (i?, Uz, 9) 
for the eigenmode shown in Fig. [3| Spectral precision is 
reached for every field. 
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Fig. 6. a)Three different conductivity profiles with 
Tbump = 1-7 (dotted blue fine), Tbump = 2.1 (solid green 
line) and Tbump = 2.8 (dashed red fine); b) corresponding 
fundamental eigenmodes (here the real part of temperature 
perturbations 9). The case with a constant radiative con- 



ductivity /vo(To) 
black line. 



-K^max is superimposed as a dot-dashed 



The first case, corresponding to the dotted (blue) line, 
denotes a hot star where the ionisation region is close 
to the surface. The conductivity hollow has an influence 
on the temperature eigenfunction which is strongly de- 
formed at its location. We will show in ii5.4lthat this de- 
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e = 0.4 and a = 1 




bump 



Fig. 7. Instability strip for the fundamental mode in the 
plane (Tbump, -4) for given values of e and a (isocontours in 
growth rates t are displayed) . The three marks correspond 
to the three particular computations done in Figs. [3] and [5] 
for Tbump = 1-7 (blue diamond, stable). Tbump = 2.1 (green 
square, unstable) and Tbump = 2.8 (red circle, stable). 



formation is not sufficient to destabilise this mode as the 
Vl'-criterion should also be taken into account (Eq. I^Tj) . 
The second one, corresponding to the sohd (green) 
line, emphasises a case where ionisation occurs approx- 
imately in the middle of the layer. The eigenvector is 
still deformed by the conductivity hollow. 
The last case, corresponding to the dashed (red) line, 
denotes a cold star where ionisation region is far from 
the surface. For this hollow location, the dynamical 
timescale is very small compared to the local thermal 
one. As a consequence, thermodynamic transformations 
are quasi-adiabatic and non-adiabatic effects are negligi- 
ble there. This is confirmed in Fig. [6)d where the eigen- 
mode is practically not deformed by the conductivity 
hollow. To enlighten this result, we have superimposed 
the case corresponding to a constant radiative conduc- 
tivity (_?i'o(To) = i^max). The obtained temperature per- 
turbations (the dot-dashed (black) line) are really next 
to the dashed (red) one, meaning that the hollow has a 
marginal effect on the mode stability. 



We then show that the conductivity hollow has only an in- 
fluence on the shape of the eigenmode in upper parts of our 
layer. In deeper regions, adiabatic thermodynamic trans- 
formations are overwhelming and the eigenvector is pretty 
unchanged by the variations in the conductivity profile. 



5.3. Parametric surveys of the instability 

As claimed in ij4.11 our conductivity profile is well-suited 
to deal with a parametric study of the K-mechanism, as 
its parameters Tbump, -A., e and a can easily be changed. 
We thus next introduce the three parametric surveys which 
allowed us to find the instability strips associated with the 
K-mechanism. 



5.3.1. The Tbump — A survey 

At first, we want to determine the influence of the hol- 
low amplitude on stability. As a consequence, we choose 
a value for a and e and make Tbump and A vary. For each 
value of these two parameters, we compute first the equilib- 
rium fields and then, the eigenvalues with their correspond- 
ing eigenvectors. Unstable modes are extracted among all 
eigenvalues as their growth rate are positive. Fig. [7] displays 
the obtained instability strip: 

— Dark (blue) areas correspond to stable fundamental 
modes (i.e. with r < 0). 

— Coloured areas correspond to unstable fundamental 
modes with the lighter the colour the bigger the growth 
rate. 

Two major results can then be derived from this fig- 
ure: (i) a particular region in our layer (Tbump G [1.8, 2.3]) 
seems to be propitious to the appearance of unstable modes, 
that is, one recovers the concept of transition region dis- 
cussed in i j3.21 (ii) a minimal amplitude in the hollow is 
required to get an instability (^min — 45%). 



5.3.2. The Tbump — cr survey 

Next, we study the influence of the slope a of the conduc- 
tivity hollow on stability. We thus choose a value for the 
amplitude A and width e while making Tbump and a vary. 
As for the previous survey, we plot in Fig.[8]the isocontours 
in growth rates t but now in the plane (Tbump, cr). 

We then found the same kind of areas than in Fig. [71 
that is, an instability strip where modes are unstable (e.g. 
the coloured region) embedded in large (dark) regions 
where modes arc stable. Nevertheless, the influence of the 
slope a on stability is less significant than the amplitude 
one. In fact, the instability strip covers the same tempera- 
ture range than in Figl?! Tbump & [1-8, 2.3], but it becomes 
almost vertical as there is no critical value of a to trigger 
off the destabilisation of the layer. In other words, there is 
a degeneracy in a as this parameter is not affecting on the 
stability. 



5.3.3. The Tbump ^ e survey 

Finally, we study the influence of the hollow width e on the 
stability by performing a survey in the (Tbump, e)-plane, 
while keeping constant A and a. Results arc displayed in 
Fig.m 

An instability strip of the same kind than the one shown 
in Fig. [7] is clearly visible as it exists a minimal value for 
the width e from which the fundamental mode becomes 
unstable (cmin — 0.15). It means that narrow hollows are 
not sufficient to initiate the instability. 

5.3.4. Summary of the surveys 

These three parametric studies allow us to show the re- 
spective influence of the hollow parameters on the layer 
stability: 

— In Figurcs[7][ni we found a particular range in Tbump, i.e. 
Tbump G [1.8, 2.3], within which the fundamental mode 
is unstable. It defines an instability strip in temperature 
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A = 70% and e = 0.4 




bump 



Fig. 8. Instability strip for the fundamental mode in the 
plane (Tbumpj cr), for given values of A and e. The blue di- 
amond, green square and red circle correspond to the same 
modes displayed in Fig. [71 
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Fig. 9. Instability strip for the fundamental mode in the 
plane (Tbump, e), for given values of A and a. 
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As the denominator is always positive, the sign of the real 
part T only depends on the numerator one. To separate the 
regions of damping from the ones of driving, we thus plot 
the real part of the work integral in Fig. [TUb . 
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Fig. 10. a) Coefficient ^ plotted over the entire box 
(Eq. I20p . The green superimposed hatched zone represents 
the location where ^ = 1 ± 0.5. b) The real part of the 
work integral numerator (Eq. [4T|l plotted for the three dif- 
ferent equilibrium models already discussed in Figs. [3] and 
ini The case with a constant radiative conductivity is su- 
perimposed as a dot-dashed black line, c) Corresponding 
radiative conductivity profiles, d) Corresponding equilib- 
rium field dK,T /dz. 



and corresponds to the particular area in the layer we 
introduced before as the transition region. 
— Only the amplitude and width of the hollow have an 
influence on the stability as we found critical values for 
both of them, i.e. Anin — 45% and Cmin — 0.15. This 
can be linked to the condition P^ which entails that 
a precise shape in the conductivity profile is needed to 
get unstable modes. Moreover, the hollow slope a is not 
a key control parameter as it does not modify the shape 
of the instability strip while varying. 

We are now going to clarify these propitious conditions 
thanks to the work integral formalism. 

5.4. Work integral in ttie non-adiabatic case 

By generalising to the non-adiabatic case the demonstra- 
tion done in Appendix I A. li we obtain the following expres- 
sion for the exact (complex) eigenvalue 



The work integral is really useful as it is possible to 
precisely check where the driving occurs. The criteria given 
in Eq. ((22l) predict that for a "sufficient" hollow, i.e. with 
a sufficient amplitude and width located in the transition 
region, the fundamental mode will be unstable. This result 
has already been checked thanks to the parametric study 
(Fig. [7] for example) . With the work integral, we consider 
the same three particular modes studied before in Figs. [3] 
and [71 

— The first one, plotted in dotted blue, expresses a case of 
a hot star where ionisation region is next to the surface. 
In this case, ionisation region is located in a place where 
density is really small, and thus \1/ ^ 1. In Fig.flOb. we 
compare this case to the one with a constant radiative 
conductivity, for which Ko{To) = i^max and dlCx/dz — 
0, and we can conclude that the conductivity hollow 
has a little influence on the work integral in this case: 
driving is unable to prevail over damping and r < 0. 
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— In the second case, plotted in solid green, the radia- 
tive conductivity begins to decrease significantly at the 
location of the transition region where ^ ~ 1 . As a con- 
sequence, driving is important in this place (Eq. [22]) • In 
addition, the radiative conductivity increase occurs in 
a place where non-adiabatic effects are already signifi- 
cant, i.e ^ < 1. This means that no damping will occur 
between the hollow position and the surface because the 
radiative flux perturbations are "frozen in" . Driving is 
overcoming damping in this case and thus, r > 0. 

— The third case, plotted in dashed red, denotes a cold star 
where ionisation region is located deeper in the stellar 
atmosphere where ^I^ 3> 1. As shown in Fig. [SI ionisa- 
tion occurs there in a quasi-adiabatic place. As a conse- 
quence, excitation provided by the conductivity hollow 
cannot balance the damping arising in the layer, thus 
r < 0. 

In conclusion. Fig. [TO] illustrates both conditions given 
in Eq. (|22|) : (i) it is first necessary to have dJCr/dz < 
to drive the oscillations. But this condition is not sufficient 
to sustain the /c-mechanism; (ii) indeed, the thermal en- 
gine underlying the conductivity hollow has to be located 
neither too deep in the star, nor too close to the surface. 
Therefore, only a specific range of effective temperatures 
allows the overlap of the transition and ionisation regions, 
which leads to the instability strips found in Figs.[7][ll 

6. Direct Numerical Simulations 

In order to confirm all of the instability strips arising in 
the previous linear stability analysis, we performed direct 
numerical simulations of the nonlinear problem. That is, 
starting from the most favourable setups found during the 
parametric surveys, we advanced in time the nonlinear hy- 
drodynamic equations to check: 

— the onset of the instability sustained by the k- 
mechanism and thus confirm the growth rates of the 
linear stability analysis. 

— the nonlinear saturation of the instability, which is of 
course not caught in the linear analysis, and have an 
estimation of the final amplitudes of modes. 

All DNS have been done using the Pencil Cod(l£|. This 
non-conservative code is a high-order centered finite dif- 
ference code (of the sixth order in space and third order 
in time) and conserved quantities are kept up to the dis- 
cretization error of the scheme. On multiprocessor comput- 
ers, it uses the MPI libraries (Message Passing Interface) 
which allow communications between processors and thus 
runs in parallel. Moreover, this code is a fully explicit one 
as the computation of the solution at time t""*"^ depends 
on the solution obtained at time t" before. The timestep 
dt is therefore limited by an usual Courant-Friedrichs-Levy 
(CFL) condition based on the consideration of the smallest 
physical timestep existing in the simulation as 



dt < min est 



Sx 



iCSty 



Sx^ 



,CStu- 



fc2 



(42) 



where est , cstx ^^nd estu are constant coefficient s depending 
on th e spatial order of the scheme (see §19.2 in [Press et al.l 
[19921 hereafter PTVF92). 



See http://www.nordita.org/software/pencil-code/ and 



iBrandenburg fc Dobieri l|2002l ) 



6.1. The needed of an implicit solver 



Equilibrium fields 




Fig. 11. Three different diffusivity profiles with Tbump = 
1.7 (dotted blue fine), Tbump = 2.1 (sofid green line) and 
Tbump = 2.8 (dashed red line). 



This CFL condition is clearly a problem in our case 
as the most favourable setups imply very large radiative 
diffusivities at the surface. Fig. [TT] emphasises this point 
where the diffusivity profile is plotted for the three com- 
mon hollows used in this work (Tbump = [1-7, 2.1, 2.8]). 
Because the ^-criterion ([^T|) imposes a weak value of Am 
and hence of the top density, one gets large surface values 
for the radiative diffusivity x = K/ pep ~ 3 in all cases. It 
means that the corresponding radiative diffusion timestep 
estx 5x'^ Ixma-x entering in Eq. P^ will be the smallest one 
and will impose a very small timestep dt. As an example, if 
we consider a typical spatial resolution of 256 x 256 (i.e. 
256 gridpoints in each direction), we obtain dt ^ 10~^ 
whereas the dynamics of the layer is rather constrained 
by the sound speed of which the dynamical timcscalc is 
fe/c^,„^^ - 3 X 10-3. 

It is therefore numerically prohibitive to reach the non- 
linear saturation of excited acoustic modes with such an 
explicit solver. To limit the number of iteration, we have 
decided to solve the diffusion of temperature implicitly. In 
fact, as implicit schemes are unconditionally stable, we have 
no more constraints on the timestep coming from radiative 
diffusion and thus the CFL becomes 



ox ox 
dt < min ( est ^cstv 



(43) 



giving dt ^ 10 ^ for the same resolution 256 x 256. 



6.2. The Alternate Direction Implicit (ADI) scheme 

In order to solve the temperat ure equation, w e adop t the 
time-split formulation given in [Malagoli et al.[ (|1995[ ) . We 
thus solve first explicitly the three hydrodynamic equations, 
i.e. density, velocity and temperature, but without solving 
the radiative diffusion term at this step. We then solve the 
temperature diffusion with the intermediate temperature 



T, 



n+l/2 



cxpl 



treated as a source term. 



The time advance of the diffusion temperature equation 
is treated implicitly in the form 



rpn-\-l rji'n 

Jt 



T, 



-1/2 



cxpl 



-T 



cxpl 



dt 



n, 



(44) 
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where the radiative diffusion term TZ is expressed by 
1 1 



7^ 



2 "+1/2 
Pcxpl *-' 



■ div [i^(T"+i)vr"+i + A'(r")VT"] 



= -[A,(T"+i) + A,(r") + A,(T"+i) + A,(r")]. 

(45) 
TZ could be directly dealt with a Cranck-Nicholson 
method and a single matrix inversion with, for instance, a 
Successive Over Relaxation (SOR) method. For a 2-D prob- 
lem, this approach forces us to invert a very sparse N'^ x N'^ 
matrix (e.g. §19.5 in PTVF92). That is why we have decided 
to implement an Alternate Direction Implicit (ADI) scheme 
which relies on operator splitting theory which is fre quently 
used i n diffusion prob lems (e.g. §19.3 in PTVF92; iDendvl 
Il977t lMasalkailll994D . It allows us to deal with tridiago- 
nal matrices (or cyclic ones depending on the boundary 
conditions used) by solving implicitly both directions suc- 
cessively. 

To treat implicitly the nonlinear terms (i.e. K{T)) 
we have adopted the following approximation com- 
monly called Rosenbroc k's method (jSaarikoski et al.lll997t 
iWitelski fc Bowenl[200l 

nn+l^l _ A /"r"~l J_ ( I (J^n+l _ rpn\ f^Q\ 



A(T"+^) = A(T") + 



J 



where J denotes the Jacobian matrix associated with the 
operator A. Expanding Eq. (j44|) in each direction thus leads 
to 



7/ \ rTin-\-l/2 rjin 

I-^jAa = A.(T«) + A,(T") + ^^eL__^^ 



I-jJ.]P 



rj.n+1 ^ T'^ + dt-P, 

. (47) 
where /, J^ and Jz are the identity matrix, a cyclic one 
and a tridiagonal one, respectively. 



6.3. Results 

All 2D-simulations were carried out using a mean resolution 
16 X 201 and a constant kinematic viscosity i/ = 5 x 10~*. 
In order to avoid as much as possible the p ropagation of 
nonradial modes (jMulet-Marquis et al.l 120071 ). we chose a 
"small" box with an aspect ratio L^/L^ = 1. Indeed, if 
we refer to classical hydrodynamic instabilities, the criti- 
cal horizontal wavelength Acrit from which the instability 
develops is generally a bit larger than the vertical extent 
Lz of the domain. As an example, Acrit/^z — '2V2 = 2.83 
in the Rayleigh-Benard convection (jChandrasekhailll96lD 
or Acrit /-^z = 2.84 in the compressible polytrope m = 1 
(|Goughet al.lll976D . 

6.3.1. Growth rates from the DNS 

Fig. [T2l emphasises the temporal evolution of the mean mo- 
mentum (puz) (where (.) is an average over the entire box) 




Fig. 12. The temporal evolution for times t £ [0, 50] of 
the mean vertical momentum (pUz) with the three different 
starting setups. In each panel, the growing (or damping) 
curve calculated thanks to the linear stability analysis is 
superimposed, a) Tbump = 2.8 (red hue), b) Tbump = 2.1 
(green line), c) Tbump = 1-7 (blue line). 

Table 1. Growth rates corresponding to the three equilib- 
rium setups given in Fig. 1121 computed in the linear stabil- 
ity analysis (the real part of eigenvalues) and in the DNS 
(computed from the exponential growth of the vertical mo- 
mentum). The corresponding relative error is also given. 



tlsb 



Tons 



Rel. err. 



red 



-2.2415 X 10" 



green 2.0484 x 10" 

blue 



2.0516 X 10"^ 1.5286 x 10" 



for the three common equilibrium setups (the red, green 
and blue ones corresponding to Tbump = [2.8, 1.7, 2.1]). 
Only one case appears to be unstable as: 

— The first (red) setup expresses the case of a cold star 
where ionisation region is deep. The linear stability 
analysis achieved before predicts that no excitation will 
occur in this case (the red circle in Fig. [7] is well outside 
the instability strip). As seen in Fig. [12k . this result 
is confirmed by the nonlinear simulation as the mean 
vertical momentum decreases with time. Moreover, the 
damping rate calculated with LSB (superimposed as a 



12 



T. Gastine and B. Dintrans: Direct numerical simulations of the K-mechanism 



red line) is reproduced with a great agreement in this 
simulation (Table [T|). 

The second (green) setup denotes the case where the 
K-mechanism is efficient. As a consequence, the funda- 
mental p-mode is expected to be unstable (the green 
square inside the instabihty strip in Fig. [?])• The DNS 
confirms this excitation (Fig. [T^b) and the growth rate 
is the same as the predicted one (Table [T]). 
The third (blue) one corresponds to a hot star where 
ionisation is next to the surface. In this case, damping 
phenomena prevail over excitation and the fundamental 
mode is stable (the blue diamond outside the instability 
strip in Fig. [7]). One more time, the result is confirmed 
by the DNS (see Fig [H; and Table [J). 



vertical profiles from the DNS and comparing them to the 
eigenvectors of the linear stability analysis. 
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Fig. 13. Upper panel: temporal evolution of the mean ver- 
tical momentum (puz). The linear stability analysis result 
from Fig.[4]is superimposed in green. Lower panel: temporal 
evolution of the corresponding envelope with the theoreti- 
cal curve of growth still superimposed. Ordinate is scaled 
logarithmically. 



In Fig. [131 we have integrated the DNS from the ini- 
tial green setup with Tbump = 2.1 till the approach to the 
nonlinear limit cycle stability, that is, the nonlinear satura- 
tion. The (green) dotted hue corresponds to the theoretical 
growth rate given in Table [J that is, t = 2.0484 x 10~^. 
The saturation of this mode appears to be achieved around 
time t ~ 200, which roughly corresponds to 170 mode peri- 
ods. Such time interval is compatible with the characteristic 
timescale of the instability given by 1/r ~ 100. Finally, the 
amplitude reached at the end of the saturation is about 
(pUz) ~ 2 X 10~^. A careful study of the modes saturation 
is beyond the scope of this paper and will be developed in 
future works. 



6.3.2. Vertical profiles 

The good agreement between the theoretical and DNS 
growth rates shown in Table [T] marks a first success in re- 
producing the K-mechanism in our simulations. We next 
address the point of the modes structure by computing the 
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Fig. 14. Upper panel: temporal power spectrum in the 
(z, u;)-plane for the radial modes k^ = 0. Lower panel: the 
resulting spectrum after an integration in depth. Dotted 
lines correspond to the theoretical frequency a;oo — 5.44 
obtained in the linear stability analysis (e.g. Fig. 2]). 



In Fig. 1141 we have performed an horizontal and tem- 
poral Fourier transforms of the vertical velocity field and 
plotted the resulting power spectrum at fc^: = in the 
plane (z, a;)-plane. With this method, we are able to 
determine exactly which acoustic modes are excited or 
not in our numerical experiment as they emerge in this 
plane as "shark fin" pea.ks ar ound given frequencies (see 
iDintrans fc Brandenbure]|2004[ ) . That is indeed what is dis- 
played in Fig. [T3] where the radial acoustic mode with a 
frequency t^oo — 5.44 is excited and the agreement with the 
linear eigenfrequency is remarkable. The second overtone 
with 0^02 — 11.05 also appears in the (z, iu;)-plane, with 
still a good agreement between the linear stability analysis 
and the DNS. In fact, one can note that this frequency is 
almost twice the fundamental mode one. It means that a 
resonance-like interaction occurs between these two modes: 
the K-mechanism excites the fundamental mode and a non- 
linear interaction transfers some energy between this mode 
and the second damped overtone, leading to the nonlinear 
saturation. 

As we have computed, thanks to LSB, the eigenfunc- 
tions from a given equilibrium setup (see Fig. [4]), one can 
compare them to the vertical velocity Uz deduced from the 
power spectra in Fig.[T31 That is what is displayed in Fig.[T51 
where a mean profile around the frequency ujqo = 5.44 has 
been computed. The DNS profile and the linear stability 
analysis one perfectly overlap meaning that we both have 
an agreement on the temporal (i.e. the frequency) and spa- 
tial (i.e. the pattern) scales of this mode. 
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Fig. 15. Comparison between normalised vertieal velocity 
profiles u^ deduced from the DNS power spectrum (solid 
line) and the linear stability analysis (dashed line). 

7. Conclusion 

We have modeled K-mechanism in Cepheids by a similar 
physical problem: the propagation of radial acoustic waves 
in a partially ionised shell. Our model consists in a perfect 
gas embedded in an entirely radiative layer and the ioni- 
sation region has been depicted by a configurable hollow 
in radiative conductivity. This approach allows to quickly 
change the shape and location of this ionised region in the 
layer through the hollow parameters Tbump, A, a and e 
(Eqs. [25ll26)) . We first have checked that a hollow with a 
sufficient amplitude and width may lead to 



dICj 
dz 



<0, 



(48) 



which corresponds to the classical instability criterion de- 
rived from the work integral and quasi-adiabatic consider- 
ations. Nevertheless, another condition is needed to obtain 
a thermal heat engine in Eddington's sense: the ionisation 
region has to be located at a certain place in the layer, nei- 
ther to deep nor to close to the surface. We have shown 
that this propitious area is located in the so-called transi- 
tion zone separating the quasi-adiabatic interior from the 
strongly non-adiabatic surface. This thermodynamical cri- 
terion can be summarised by 



T; 



lonisation 



Tr 



TR 



and 



(ci,TTR,)Am 



1. 



(49) 



With this second criterion, appropriate boundary con- 
ditions have been chosen in order to have a transition re- 
gion in the middle of our box for the fundamental acoustic 
mode. Both radiative and hydrostatic equilibria have then 
been discretized on the (spectral) Gauss-Lobatto grid and 
solutions have been computed thanks to a linear solver. 
Once done, we have performed a linear stability analysis by 
computing the whole spectrum as well as associated eigen- 
vectors for the radial oscillations equations. 

The main advantage of our approach lies in its ability to 
make the parameters of the conductivity hollow vary. The 
major result of parametric surveys is the checking of the 
previous conditions by way of the appearing of instability 
strips: only configurations satisfying both conditions (j48t - 
HU]) have led to unstable fundamental modes (e.g. Figs. [7] 
and[TU|. 

Then, these most favourable setups found in the linear 
stability analysis have been the starting points of our first 
2-D (nonlinear) DNS of K-mechanism. However, because of 
a too restrictive Courant-Friedrichs-Levy (CFL) condition 



at the surface, these setups were not well-suited for a fully 
explicit code such as the Pencil Code one. We therefore 
have developed a new implicit module to deal with the 
radiative diffusion term and thus soften the CFL constraint. 
These DNS have confirmed with a great agreement both 
growth rates and structures of linearly unstable modes (see 
Figs. [T51 and [T51) . Moreover, we have been able to reach the 
nonlinear saturation which involves an intricate coupling 
between the fundamental mode and (at least) the damped 
second overtone of which the period is a multiple of the 
fundamental one. 

This work constitutes a first step in our Cepheids' 
project devoted to the modelling of the convection- 
pulsation interaction in the coldest Cepheids close to the 
red edge of the instability strip. Indeed, it is well known 
that convection occurs in a non-negligible part of these stars 
and modify the pulsation properties. For now, models based 
on time-dependent convection theories reproduce quite well 
the position of this red e dge or, say, obser vational periods 
and hght curves (YKB98: lBono et al.|[T999l ). However, they 
involve many free parameters which are either fitted to the 
observations or hardly constrained by theoretical values, 
leading to almost similar results for diff erent parameter sets 
(iKoUath et al.ll2002l : ISzabo et al.ll2007D . 

Despite their own limitations (their weak contrasts in 
pressure through the computa tional domain or their the r- 
mal timcscalc problem, see e.g. lBrandenburg et al.l (|2000( )). 
DNS are without a doubt a good way to address this 
convection-pulsation interaction as they fully take into ac- 
count the crucial nonlinearities. With a conductivity pro- 
file solely based on temperature (Eq. E5)) . it is easy to lo- 
cally shape one (or several) convective zone by increasing 
the temperature gradient above the adiabatic one. Indeed, 
as an ionisation region corresponds to a local increase in 
opacity, it is well known that convection can develop there. 
This occurs in cold Cepheids where two separate convec- 
tively unstable regions superimpose with the Hel/H and 
Hell ionisation regions. By adjusting in DNS the convec- 
tive zone width or the strength of gravity, it will be possible 
to match the (local) turnover timescale of convection with 
the mean mode period, and then to study the coupling be- 
tween convection and pulsation. 
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Appendix A: Work Integral 

A.l. Using the quasi-adiabatic approximation 

In the energy equation in Syst. ([3]), we can first separate 
adiabatic terms from non-adiabatic ones as 






1 



■divF', 



POCyX 

and, by using the ideal gas equation of state H]), 

7-1 



P — Padia 



A 



■divF'. 



(A.l) 



(A.2) 



This equation corresponds to the pressure perturbation due 
to both adiabatic and non-adiabatic oscillations. The mo- 
mentum equation in Syst. ([3]) then becomes 
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{X + SX){u + Su) = Vp'^^., 

Po 



Padi', 

Po 



Adiabatic part 

1 7- 1 



(A.3) 



Po A 



V(divF') 



Non-adiabatic perturbation 



This equation can formally be written as a generalised 
eigenvalue problem with a perturbation operator 



(A + 6A){iP + 5xP) = (A + SX)B{'ijj + (5-0), 



(A.4) 



where SX and S^p are caused by the perturbation operator 
6 A (here the non-adiabatic effects). 5X can then be obtained 
from a first-order perturbation an a,lysis equivalent to these 
done in quantum mec hanics (e.g. iBender fc Orszaeill978t 
iDvson fc Schutzlll979h 



SX^ 



(ipolBtPo) 



(A.5) 



where the symbol ( | ) means t he following dot product (e.g. 
iLvnden-Beil fc Qstrikeilll967f) 



(/i|/2)= / fl-f2PodV. 
Jv 



(A.6) 



We thus obtain the expression for the eigenvalue perturba- 
tion SX 



6X 



[ M*.V(divF')dV^ 
1 Jv 



\u\^PodV 



(A.7) 



We can go a step further using 



M*.V(divF') = div(M*divF')-divM*divF', (A.8) 

of which the first term in the RHS vanishes due to the rigid 
boundaries for the velocity as 

f div {u* div F') dV = ^u* div F' -(13 = 0. (A.9) 

The second term in the RHS can also be simplified by using 
the continuity equation 



divM* = —A 



JP* 



PO 



(A.IO) 



where Sp is the Lagrangian perturbation of density. We fi- 
nally get the expression for the eigenvalue perturbation in 
the non-adiabatic case as 



SX = - 



[ {j-l)^dWF'dV 
Jv P 



H^PodV 



(A.ll) 



where we have assumed that the adiabatic eigenmode is 
purely imaginary, that is A = iw, and thus A*/A = —1. We 
note that 5ft((5A) corresponds to the damping (or growing) 
rate of a non-adiabatic mode, commonly written r. 



A. 2. A thermodynamical approach 



It is also possible to derive Eq. (jA.lip by energetic con- 
siderations. Her eafter, we will essentially follow Hansen's 
demonstration ([Hansen fc Kawa ler 199j). Let us adopt a 
thin shell of mass dm at a certain radius. During a com- 
plete cycle of oscillations the work dW done by the shell 
on its surroundings is linked to the internal energy dU and 
heat dq gained by this shell through the first principle of 
thermodynamics 

dq^dU + dW. (A.12) 

Integrating over a cycle of oscillations leads to 



dq 



dU 



dW. 



(A.13) 



We are now going to suppose that the oscillation mecha- 
nism is quasi-adiabatic. It implies that every thin shell dm 
behaves as a Carnot-like heat engine where each process is 
reversible and the shell comes back to its initial position 
after each cycle of oscillations. Internal energy dU being a 
state variable, one gets § dU ~ hence 



W ^ idq. 



(A.14) 



Now one applies the second principle of thermodynamics 
claiming that 



dS 



dq 



(A.15) 



After some time elapsed, we have T ~ Tq + ^T and then 

(A.16) 



db ^ — - -2 dq. 



as a first-order approximation, with § dS = for the same 
reason than U . We thus write 



dq= <t> jrdq. 



(A.17) 



With Eqs. (jA.14[) and (|A.17[) . the work over a cycle is given 
by 



W=<l> ^^dq, 



(A.18) 



which leads to the total work integrated over mass shells 

6T 



Wt, 



M 



Tr 



■dqdm. 



(A.19) 



If VFtot > 0, the star produces work over a cycle of oscil- 
lations and the initial perturbation will grow. In this case, 
the star is driving pulsation which is thus unstable. It may 
occur when shells gain heat (i.e. dq > 0) during compres- 
sion phases (i.e. 5T > 0) and this is the s o-called "valve- 
mechanism" proposed bv lEddingtonl ()1917f ). 
Let us write the energy equation in the form 



dq = divF'di, 

Po 



(A.20) 



and substitute it in Eq. (jA.19[) 
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dWtc 
dt 



XT' 

— divF'dV. 

V ^0 



(A.21) 



Kinetic energy over a period of oscillations is given by 

E=^J \u\'p„dV, (A.22) 

and Eqs. (|A.2mA.22|) lead to the following expression for 
the growth rate of a mode 



1 dWtot/dt 



(A.23) 



Note that the factor 1/2 comes from the fact that the 
pulsation amplitude grows or decays one half as docs en- 
ergy. The two formulations Eqs. (jA.lip and (|A.23p are of 
course alike by assuming the quasi-adiabatic approximation 
ST/To-{j-l)Sp/po. 

Appendix B: The radiative term in energy equation 

Radiative flux perturbations can be written as 

- div F' = div [KqVT' + K'VTo) . (B.l) 

By assuming 9 = T'/Tq and using Eq. ^, it entails 



d , f dO dTn dTf) , 

— F' = div KoTo— + Ko9-^ + K^KtO-^ ) e, 
dz V dz dz dz 



We then expand /Ct with equilibrium equations (|29p as 
d In A'o To dA'o 



ICi 



and thus obtain 



- div F' = — 
dz 

Finally, one gets 



divF' 



din To Fbot dz 



(B.2) 
as 

(B.3) 



A'oTo^ + i^o^^+^To^). (B.4) 
dz dz dz 



^ [K^T^e) = A, (XoTofl) . (B.5) 
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